Quantum simulations of optical systems 
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Within a framework of a two-dimensional microscopic purely-quantum mechanical model we 
analyze dynamics of single-photon wave packets interacting with optical elements (beam splitters, 
mirrors) modeled as systems of two-level atoms. That is, we utilize a two dimensional cavity to 
simulate the quantum behavior of simple optical components and networks thereof. The field is 
quantized using the canonical procedure, and only the basis states with one unit of excitation 
are included. This, however, covers the linear optical phenomena. The field is taken to interact 
with localized atoms through a dipole interaction. Using different configurations of atoms and 
choosing their frequencies to be resonant or off-resonance, we can model mirrors, beam splitters, 
focusing devices and multicomponent systems. Thus we can model arbitrary linear networks of 
optical components. We show the time evolution of a photon wave packet in an interferometer as an 
example. As the state of the field is known at each instant, spectral properties and spatial coherence 
can immediately be obtained from the simulations. We also know the states of the two level atoms 
constituting the components, which allows us to consider their quantum behavior. Here the decay 
of an excited atom into the vacuum state of the electromagnetic field in the two-dimensional cavity 
is studied. 

42.50.Ct,32.80.-t,32.90.-l-a 



I. INTRODUCTION 

It is well understood that the electromagnetic fields 
giving rise to all optical phenomena have ultimately to 
be represented by quantum operators. These couple to 
the degrees of freedom of matter and their modification 
due to this interaction constitutes the quantum counter- 
part of the action of optical components. Ordinary op- 
tical devices operate in the linear regime of interaction, 
but the important area of Nonlinear Optics is based on 
higher order effects of the field-matter coupling. In most 
situations, the optical phenomena can be described en- 
tirely in terms of classical fields, but many recent in- 
vestigations require that the quantum character of the 
field is accounted for. Such research constitutes the top- 
ics of Quantum Optics [Q-g|. However, many quantum 
effects are of interest even in the linear regime of opera- 
tion: quantum noise [H sets the limit to communication 
by optical channels and amplifiers, quantum interference 
shows up in precision measurements and tests of funda- 
mental issues and also in reading and writing quantum 
information. Manipulation of quantum information such 
as quantum computations usually requires the inclusion 
of higher order effects, i.e. nonlinear interactions between 
the qubits ^. 

In this paper, we are going to discuss the dynamics 
of single photon wave packets in various two-dimensional 
atomic configurations. These are taken to be models of 
optical networks, where we explicitly include the atomic 
nature of the optical components distributed over the 
volume under investigation. This approach provide us 



with a completely microscopic quantum-mechanical pic- 
ture of how photon wave packets interact with optical 
elements represented as collections of two level atoms. 
For practical reasons, we have to restrict our work to 
one-photon states, but this is not such a serious limita- 
tion as it may seem. All linear optical effects are based 
on the single-photon interacting with material structures, 
and consequently we have a general description of Quan- 
tum Optics phenomena in the linear regime. The need 
to consider multi-photon effects arises only in connection 
with the quantum treatment of Nonlinear Optics. 

There are two basic ways to approach the quantization 
of optical systems. In the conventional one, we deter- 
mine a complete set of eigenmodes of the total universe, 
and express the fields of interest in terms of these. Any 
matter present is described through its interaction with 
the fields, and the coupled field-matter problem is then 
solved to the best of our ability. This is the approach 
utilized in traditional QED, and its development is found 
in many standard texts. The alternative approach, de- 
signed for Quantum Optics applications, is to determine 
the eigenmodes of the system at the classical level, the 
matter involved is then treated as boundary conditions 
on the field modes. Especially the new area of Cavity 
QED research ra] , utilizes this point of view, and it pro- 
vides the basis both for quantum communication theories 
and many fundamental investigations. 

In the field of optics, the components are usually 
treated as boundary conditions only, and the complete 
optical device is considered to be an optical network. 
This approach has been discussed thoroughly in the clas- 



sical regime of operation 0. For linear devices, the 
classical treatment can be taken over into the quan- 
tum regime by the use of suitable Quantum Optics tools 
P-[lO|. In principle, any device understood classically, 
can be treated quantum mechanically with such an ap- 
proach. The specific quantum features manifest them- 
selves in the initial conditions and the restrictions on ob- 
servability imposed by quantum theory |ll| . 

Another specifically quantum mechanical effect is 
the occurrence of spontaneous decay. Within a one- 
dimensional model of the modes of the universe, this is 
discussed in Ref . |]T2[ , where both free Weisskopf- Wigner 
decay and cavity modified decay are discussed. Such phe- 
nomena have been the object of much interest within 
tED research, for an extensive list of references see Ref. 
2|. Within the model chosen there, one can see the 
emergence of the exponential law and the inhibition of 
decay observed in a photonic band gap structure. In gen- 
eral, the model provides insight into the role of atomic 
media in the irreversible transfer of excitation energy into 
the field modes of the universe. 

In this paper we combine the two views discussed 
above: We retain a description in terms of a complete set 
of two-dimensional eigenmodes of the universe. The op- 
tical components are described in terms of their atomic 
constituents. All atomic structures are represented by 
spatially localized two-level atoms. These are treated as 
point-like partciles in accordance with the dipole approxi- 
mation assumed to be valid. The state of the field is taken 
to be a single photon wave packet with a narrow energy 
distribution. In this case, the state can be described by 
a truncated expansion in terms of modes of the universe. 
The spatially distributed two-level atoms describing the 
structures are taken to be initially in their ground states. 
The atoms can be chosen to resonate with the central 
frequency of the photon wave packet or be well off reso- 
nance; various effects can be modeled in this way. When 
the single photon is absorbed, only one of the atoms is ex- 
cited, and the field is reduced to its ground state. Such a 
choice limits the Hilbert space needed in the calculations 
to manageable size, but allows us to investigate many 
simple networks of significance in linear optics. All such 
effects are, in principle, describable at the single photon 
level; only Nonlinear Optics effects require more photons, 
which would make the Hilbert space expand beyond the 
limits of available computer resources. 

Our approach based on a complete set of eigenmodes 
allows us to investigate the dynamic performance of many 
linear systems. In order to illustrate the method, we 
select the simplest optical components: mirrors, beam 
splitters, focusing devices and interferometers. The over- 
all performance of the components follows directly from 
their classical theory, but our approach allows us to inves- 
tigate the microscopic (quantum) behavior of the setup. 
Quantum coherence between various spatial regions in 
the device is directly visible in the states calculated, and 
the time and space scales of the various interferometric 
structures can be read off the results. Combined with 



various models of measurements, our calculations con- 
tain considerably more information than a simple classi- 
cal computation. Here we only discuss the measurement 
of frequency and the possible occurrence of filtering ac- 
tion in the atomic structures, which does not in itself de- 
pend too much on the quantum nature of the fields. But 
modeling the frequency detection by atomic absorption, 
we utilize the full character of the model, which allows 
further extension to quantum correlation measurements 
if we so desire. 

Our work is based on the model put forward in 
|13| which we extend to two dimensions. The quan- 
tized modes of the universe are introduced in Sec. II 
together with their interaction with the spatially dis- 
tributed atoms. In Sec. Ill we specify the details of 
the model and indicate how the calculations have been 
carried out. Section IV presents the various simple com- 
ponents analyzed in this paper. We describe how they 
are modeled and show the results of the detailed solution 
of the time evolution. Finally in Sec. V we present our 
conclusions and discuss possible extensions and applica- 
tions of the work. 



II. OPERATORS FOR THE FREE FIELD IN TWO 
DIMENSIONS 

The field is enclosed inside a two dimensional cavity 
determined by the relations 
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(2.1) 



The periodic boundary conditions restrict the allowed 
values in k-space to a discrete set 
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(2.2) 



In computer simulations, the k-values must be restricted 
by giving some upper limit for the integer n which cor- 
responds to a specific frequency cut-off. The electric and 
magnetic field can be expanded [|| using the mode func- 
tions 
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B(r) = } E f ^^ ) («k.(k X ek.)e^''-- - h.c), (2.4) 



where the summation ^ is over all k-values (2.2) and 

ks 

two polarization indices s — 1,2. The frequency Wks is 
the same for both polarizations 

Wks = c|k|. (2.5) 

The general k-vector in two dimensions can be written 



k = k^ei + kye2 = |k|(cos((/))ei + sm((/>)e2). (2.6) 

The polarization vectors which obey the usual right hand 
rule conventions are 
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k X Eki = -kyCi + kxe2 
k X ek2 = |k|e3. 



ek2 = - sin(0)ei + cos{4>)e2- 



(2.7) 

f2 g\ The energy-density operator is 



The k-vector and polarization indexes satisfy the rela- Hir) — —e E^fr) + — — B^fr) 

tions iP 2 2^0 
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(2.9) Using (2.3) and (2.4) gives 
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In our simulations, we have restricted the polarization 
of the field to eki . The modes with s = 2 are taken to 
have have zero amplitudes. In addition to that we restrict 
the number of excitations of our basis vectors to one. For 
these kind of basis vectors the terms aksSk's' and a\^gO]^isi 
do not give any contribution. These terms can be omit- 
ted from the expressions. For the states described above, 
the expectation values are obtained by replacing the op- 
erators with the coefficients of the corresponding stat- 
evectors ak -^ Ck and aj^ -^ cj^. The normally-ordered 
terms in the energy density become (normal ordering is 
indicated by colons) 
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The two-fold summation over the k-space is seen to fac- 
torize and the formulas for R and 5^ are Fourier trans- 
forms of two different functions. For numerical simula- 
tions these two properties are essential as will be seen 
later. We note that if the polarization is such that the 
modes with s = 1 are taken to have zero amplitudes, 
then the two terms : ieoE2(r) : and : ■^&{r) : in the 
expression for the energy density are equal. 



Integrating ( |2.13 ) over the spatial coordinates and us- 
ing the integral 
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gives the familiar form 
i^F = - ^ fiuJk{dl_ak + flkftk) = ^ huJi^{aldk + -), 
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(2.21) 



which in the normally ordered form reads : Hp 
J2 Tit^kakSk 



III. THE GENERAL HAMILTONIAN AND THE 

STATES 

In the previous chapter the formulas for the field in the 
vacuum were derived. In this chapter we add an assembly 
of two level atoms to the cavity and give the correspond- 
ing Hamiltonians. The general form of the statevector 
with one excitation is also given. The material presented 
here is based on the similar simulations in one dimension 
done by V.Buzek et.al. ||lj]. The simulations in two di- 
mensions are numerically more demanding, but we have 
been able to develop efficient numerical methods which 
make these simulations possible. 



A. The Hamiltonian 

The total Hamiltonian H can be divided into three 
parts 



H = Hf + Ha + Hi, 



(3.1) 



where the field Hamiltonian is given by equation ( 2.21 ). 
The atomic Hamiltonian is the sum over all one-atom 
Hamiltonians 
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(3.2) 



where ldj is the transition frequency of the j-th atom 
and (7^ is Pauli's spin matrix. In the interaction Hamil- 
tonian the dipole approximation is used. For simplicity 
the dipole operator is taken to be 



T>, = {D,a\+D*al)h. 



(3.3) 



i.e. it has a component in the 63 direction only. The 
general dipole vector would have components in the x- 
and y-directions too. The interaction Hamiltonian has 
the form 



Na / Na 

Na 

k j=i 

The first sum contains all the basis vectors where the 
excitation is in one of the field modes and all the atoms 
are in the ground state. In the second sum the field 
modes are in the vacuum state and one of the atoms is 
excited. The complex numbers Ck and Cj are the proba- 
bility amplitudes of the corresponding basis vectors. We 
have dropped the polarization indices because in our sim- 
ulations only the basis vectors with the polarization vec- 
tor Eki arc excited as was discussed earlier. 

The general Gaussian one photon statevector is of the 
form 



i*) = E^kiik,{o}), 



(3.8) 
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Hi = 






(3.4) 



where E(rj) is the electric field operator ( |2.3D at the po- 
sition of the atom. The rotating wave approximation 
(RWA) is to be used, and we neglect the a+^a'' - and aia 
terms. In addition to that we replace the mode frequency 
in the electric field operator by the atomic frequency and 
use the dot products 63 • Cki = — 1 and 63 • ek2 = to get 
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(3.5) 



in what follows we omit the polarization index in sub- 
scripts of field operators. The coupling constant is 
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(3.6) 



Only those modes whose resonance frequency is close to 
the atomic frequency interact significantly with the atom, 
so we can replace the mode frequency Wk by the atomic 
frequency luj in equation (3.6). 



B. The statevector 

In all simulations we have restricted the total number 
of excitations to one. Consequently, the most general 
statevector of the atom-field system has the form 



l*> 




where the mode coefficient Ck is 

Ck = , exp ( TTrikx — k^n) 



(3.9) 
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The parameters M and A? , are 
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If the cross-variance A^^ , vanishes the formula for Ck 
reduces to two independent Gaussian distributions 



Ck = (2^AL)~^/'(2^A2 
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All ini tial d istributions used in o ur si mulations are of the 
form ( |3.12| ). The distribution ( 3.12 ) in k-space is cen- 



tered around {k^o, kyo) with the corresponding central 
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A^ the distribution is symmet- 



frequency ojQ' 

ric. If A^.^ < A^ the distribution is wider in the y- 
direction (and vice versa). The variances in k-space and 
configuration space are inversely proportional. If Af.^ 
is small, the energy density distribution in configuration 
space is wide in the x-direction. The normal ly o rder ed en - 
ergy distribution associated with the state (3.£) or ( 3.12| ) 
is well localized near the point ro in the configuration 
space. Essential for this is the phase part e~^^''^° of the 
coefficient Ck. If the form of the phase was different the 
intensity profile would not be Gaussian. 

The time evolution of the Gaussian wave packet in- 
side an em pty cavity is determined by the Hamiltonian 
Hp (2.21) with the corresponding evolution operator 



exp{—j-HFt). Applying this to the state ( |3.7| ) gives for 
the time-evolution of the coefficients Ck(t) = Ck(0)e~*"''*. 
The absolute value of the coefficients remain the same, 
only the phase changes. For the phase part we get 

exp(— ik • Tq — iuji^t) = exp[— ik • (ro + cte^)], (3.13) 

where k = |k|ek. The time evolution inside the empty 
cavity reduces to the time evolution of the parameter 
r{t) = To + cteic. We remember that the phase factor 
determine the shape of the normal ordered intensity pro- 
file. Because the time evolution of the phase is different 
for different modes, the normal ordered intensity does 
not preserve its original Gaussian shape. If the direction 
of the vector et is more or less the same for all basis 
vectors which have nonzero coefficients, the shape of the 
energy density distribution remains approximately the 
same longer. The situation is like this when the state 
vector in k-space is centered around some k-value far 
from the origin and the variances are small. 



D. Numerical methods 



1. Integration of Schrodinger equation 



Our choice for the integration method of the time de- 
pendent Schrodinger equation is a classical four stage 
fourth order Runge-Kutta method. If the wavefunction 
at time t is \'^{t)) the wavefunction at a later time t + At 
(At small) |^(i)) is given by the following algorithm |1J] 



|fci) = AtH\^it)) 
\k2)^AtH{\^it))+0.5\ki)) 
\k3)^AtH{\^{t))+Q.5\k2)) 
\k^)=AtH{\^{t))+0.5\h)) 



|*(t + At))^|*(t)) + M + M + M 



(3.19) 



^ + 0((At)^) 



C. Transformation to the interaction picture 

It turned out to be faster to carry out the numerical 
integration in the interaction picture. The transforma- 
tion Hamiltonian is Hq = Ha + Hp. The interaction 
Hamiltonian in the rotating frame is 

ijf ' = cxp{iHot/h)Hi cxpi-iHot/h), (3.14) 

which is obtained by the following replacement 
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in equation (|3.5D, and we get 
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The statevectors in the interaction picture become 

|1')(^) = exp(iiJoV?i)l*> (3-17) 
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Ecke-''*|lk,{0})+5]c,e-^*|{0},l, 

k j = l 



and the Schrodinger equation for the wavefunction is 

(3.18) 



ih- 



dt 
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Integration of the Schrodinger equation in the interaction 
picture is faster than the original equation because only 
the interaction Hamiltonian is present. 



The timestep At is a fixed constant. 

The essential part of the integration from the numeri- 
cal point of view is how to evaluate the right hand part 
of the equation (3.18) as eff iciently as possible. The first 
term in the equation (3.16) gives 
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and the second one 
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Hence the new coefficients for the atomic (c' ) and field 
(cjj.) basis vectors become 
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where 
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Both T{r,t) and C/(k, t) are two dimensional Fourier 
transforms, so in numerical calculations the Fast Fourier 
Transform (FFT) can be used. The speed increase ob- 
tained by using FFT instead of the direct summation is 
enormous especially in simulations with a large number 
of atoms. In some simulations it can be said that only 
this method makes these simulations possible. 

There are several natural checks for the numerical sim- 
ulations. First of all the norm of the wavefunction has 
to remain unity for all times. The system is closed so 
the total energy of the system must be constant all the 
time. Th e fiel d energy can be calculated using either the 
formula ( 2.21 ) or integrating the energy density over the 
whole cavity. The two methods should give the same 
results. 



2. A method to detect a local time dependent spectrum 

In the following simulations the spectrum is detected 
using so-called analyzer atoms [Q. Many atoms with 
a very small dipolc coupling constant are put into spe- 
cific locations in the cavity. All the atoms have different 
transition frequencies in between LOmm and Umax 
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j = l,2...N 



(3.28) 



Also the dipole constants are all different and very small 

C 



Dj = — , 



(3.29) 



where C is a ve ry sm all constant, typically C=0.0001 or 
so. The form (3.29) of Dj gives the same decay con- 
stant r for all the atoms because in two dimensions F 
is directly proportional to the product D'juj'j. Because 
the dipole coupling is small, the atoms have very small 
decay constants and linewidths and only the radiation 
which is exactly on resonance with the atom can excite 
it. Therefore the excitation of the atoms as function of 
Lo can be interpreted as a spectrum of the field at the 
position of the atoms. Because the interaction between 
the radiation and the atoms is small, the state of the field 
does not change appreciably. The method can be used 
to detect the local time dependent spectrum. Two-time 
averages, usually used in spectrum calculations, are not 
needed. A more detailed description of the method and 



comparisons with the time dependent spectra defined us- 
ing two-time averages |16| can be found in the paper by 
M.Havukainen and S.Stenholm |Iq], where it was used 
to detect the spectrum of a radiation emitted by a laser 
driven three level atom. 



IV. SIMULATIONS 

In this chapter the results of several simulations are 
presented. First we show that the energy density profile 
of the free photon does not preserve its shape if ojq is 
small as was explained earlier. In the second and third 
simulation, atoms are used as mirrors and beam splitters. 
Using these components it is possible to build many opti- 
cal systems. We present an interferometer as an example. 
We also present a simulation of a two-slit experiment. Fi- 
naly, we also briefly study a spontaneous decay of a two 
level atom into the vacuum of electromagnetic modes in 
the two-dimensional cavity. 



A. A free photon 

In the first simulation the time evolution of the free 
photon wave packet is studied. The initial wave packet 
is Gaussian (3.12) with parameters xq = —8.0, ya — 0.0, 



4.0, h 



yo 



0.0 and A|^ 



^ly 



1.0. The prob- 



abilities |ckp of the field modes are shown Fig. ^. The 
central frequency of the photon wave packet is so small 
that the k- vectors of the modes with nonzero amplitudes 
are not parallel. We would expect this to be observed 
as was explained earlier. The time evolution of the en- 
ergy density at two time values is shown in Fig. g. The 
initial Gaussian photon wave packet has an energy den- 
sity centered at x — —8.0, y — 0.0. The wave packet 
is moving to the right. During the free evolution energy 
density becomes delocalized. From the figure we see that 
at i = 20.0 the width in the y-direction is much larger 
than the initial value. This spread of the width of the 
original wave packet is a standard quantum-mechanical 
effect. 



Initial state 



B. A mirror 



0.04 



0.02 




FIG. 1. The probabilities |ckp of the Gaussi an in itial 
state in k-space. The parameters in the equation (pT3) are 
xo = -8.0, yo = 0.0, fc^o = 4.0, kyo = 0.0, A^ = 1.0 and 
A|y = 1.0. Here we consider the size of the cavity to be 
L = IOtt and we take into account 256 x 256 modes of the 
electromagnetic field. Only one polarization (s — 1) is taken 
into account. 



It is possible to "build" mirrors and beam splitters us- 
ing two level atoms. In the next simulation many atoms 
with large dipole constants were arranged into a slab con- 
figuration. We take a 45° angle between the slab and the 
X-axis. We assume all atoms to have the same transi- 
tion frequencies and dipole constants. Th e init ial photon 
wave packet has a Gaussian distribution ( 3.12 ) with pa- 
rameters kxo = 5.0, kyQ — 0.0 and A|^ — A'ly — 0.125. 
The atoms in the slab are exactly on resonance with the 
incoming photon wave packet (i.e. the central frequency 
of the wave packet luq = 5.0 is equal to the transition 
frequency of the atoms). The dipole constant is large 
D — 0.5. We assume that the mirror is composed of eight 
layers of atoms as close to each other as possible. In our 
case we assume to distance between neighboring layers 
of the atoms to coincide with the grid in the configura- 
tion space (the grid spacing is Aa; and for the given ori- 
entation of the mirror the distance between the different 
atomic layers is chosen to be AX = a/2- Ax = 0.17). The 
central wavelength of the incoming photon wave packet is 
A — 1.26 so the difference between the neighboring atoms 
much shorter than the wavelength of the incoming wave 
packet. 

We plot the energy density of the one-photon wave 
packet reflected by the mirror in Fig. 0. Firstly we plot 
the initial wave packet at t — 0.0 (ajT The photon is 
coming towards the atoms of the mirror. These atoms 
become excited by the incoming wave packet. The "sec- 
ondary" radiation which is emitted by the atoms interfere 
with the incoming wave packet. This secondary radiation 
can formally be expressed as a sum of the two terms - 
the first destructively interfere with the incoming wave 
packet. As a consequence of this interference the incom- 
ing wave packet is "destroyed" (i.e. becomes extinct). 
The other part of the radiation which is "collectively" 
radiated by the atoms of the mirror represents the re- 
flected wave packet. In fact, the process of reflection of 
the wave packet by atoms of the mirror represents purely 
quantum (microscopic) version of the Ewald-Oseen ex- 
tinction theorem [17|. In Fig. g(b) we have chosen con- 
ditions such that at t = 20.0 all the radiation is reflected 
by the atoms. 



FIG. 2. The time evolution of the energy density of of the 
initial Gaussian photon in free space. The parameters are the 
same as in Fig. hi We see that the initial wave packet (a) is 
nicely localized in the configuration space while at later times 
it does not preserve its initial shape - we see (b) the spreading 
of the original wave packet in the t/-direction. 



C. A beam splitter 
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FIG. 3. The energy density of tlie one-plioton wave packet 
reflected by a mirror co mpose d of two-level atoms. The ini- 
tial photon is Gaussian (3.12) with parameters a;o = —8.0, 
yo = 0.0, fc^o = 5.0, kyQ = 0.0, AL = 0.125 and A|y = 0.125. 
The atoms of the mirror are exactly on resonance with the 
central frequency of the photon wave packet [lo = 5.0). The 
dipole constant of the atom is chosen to be D = 0.5. The total 
number of atoms considered in this simulation was 1584. The 
number of modes is the same as in the simulation presented 
in Fig. |. 

The direction of propagation of the reflected wave 
packet is the same as expected in the classical the- 
ory. The energy density compared to the incoming wave 
packet is changed but is still clearly localized. Note that 
the energy density is not perfectly symmetrical. The rea- 
son is the same as in the simulation with a free photon, 
i.e. the distribution in fc-space is broad and near the 
origin so the spread of the wave packet is clearly seen. 
Additionally, the interference between components of ra- 
diation emitted by different atoms of the mirror plays a 
role. In the left part of Fig. ^ we see the energy den- 
sity of the photon wave packet close to the surface of the 
mirror. We see that the incoming and reflected parts in- 
terfere. We also see that no energy is transmitted by the 
atomic slab. In this sense the atoms serve as a mirror. 
Nevertheless, one has to remember that the atoms dur- 
ing the process of reflection of the original wave packet 
become excited, that is the mirror under consideration 
has its own "internal" (quantum) degrees of freedom, so 
the part of the original energy can be (transiently) ab- 
sorbed by the mirror. This also result in the fact that 
this quantum mirror might become entangled with the 
reflected wave packet. 

We note that the parameters of the atoms in this simu- 
lation were carefully chosen in such a way that the atoms 
really form a mirror. If the parameters are changed then 
part of the radiation can be transmitted, that is the col- 
lection of the atoms can play a role of a beam splitter. 



In the previous simulation we have shown that it is 
possible to build an almost perfect mirror using two level 
atoms, assuming the parameters are chosen correctly. Us- 
ing slightly different parameters, we flnd that the atoms 
can behave as a beam splitter. There are several ways 
how to modify the "mirror" configuration to obtain a 
beam splitter - for instance, we can consider a smaller 
number of atoms, or we can decrease the dipole con- 
stants, or change the resonance frequencies of the atoms. 
We tried all the possibilities and the most satisfactory 
results were obtained by detuning the atoms. The fre- 
quencies of the atoms are now taken to be w = 10.4. The 
center frequency of the incoming photon wave packet is 
loq = 15.0, i.e. the detuning is really large. The time evo- 
lution of the energy density of the electromagnetic field 
in this case is shown in Fig. 0. The line in the middle 
represents the positions of the detuned atoms. There is 
only one layer of atoms instead of eight as in the mirror 
simulation. At t = 0.0 the photon is propagating towards 
the atoms. Here again the incoming wave packet excites 
the atoms. Now the quantum interference between the 
incoming and emitted radiation is such that part of the 
original wave packet is transmitted by the layer of the 
atoms. The other part is reflected. In Fig. ^ we clearly 
see that at i = 20.0 the original wave packet is split into 
two parts propagating up and to the right. The energy 
is divided equally, that is the atoms form a 50-50 beam 
splitter for the incoming photon. Here we stress that the 
beam splitter under consideration has its own internal 
degrees of freedom and transiently it becomes excited. 
Nevertheless, after a while the atoms completely emit 
the excitation energy and the beam splitter is in a ground 
state - at this point it is completely disentangled from the 
one-photon radiation held which is now in a pure super- 
position state with two macroscopically distinguishable 
components (reflected and transmitted). 



FIG. 4. The energy density of the photon wave packet 
which is sphtted by a quantum beam sphtter composed of 
a set of two level atoms composing a one-dimensional crystal 
- the quantum beam splitter. The initial photon is Gaussian 
(3.12) with parameters xq = —10.0, j/o ~ 0.0, k^o = 15.0, 
kyo = 0.0, AL = 0.125 and A^^ = 0.125. The transition 
frequency u — 10.4 of the atoms is detuned from the central 
frequency of the incoming photon wave packet uq = 15.0. The 
total number of atoms is equal to 881, while the number of 
modes is the same as in the simulation presented in Fig. bl 
Here we again assume the dipole constant of the atoms to be 
D = 0.5. 



Photon at the mirror 



Photon at the beamsplitter 





In the right hand part of Fig. |^ energy density of the 
photon wave packet is shown close to the "surface" of 
the beam splitter. To the left of the atoms the incoming 
and reflected wave packets interfere. We also see that 
a fraction of the original radiation is able to "pass" the 
atoms and to continue to propagate to the right. The 
wavelength of the photon was chosen shorter than in the 
mirror simulation, which can be seen from the interfer- 
ence structure. 

We have also studied spectral properties of reflected 
and transmitted parts of the original wave packet. In 
this situation 200 atoms were used to detect the time de- 
pendent spectra of the two outgoing parts of the photon 
by applying the method described earlier. Both spectra 
were identical to the spectrum of the incoming photon. 
This means that our quantum beam splitters and mirrors 
are linear devices, which is important if we want to build 
optical networks out of the considered optical elements. 



D. Parabolic mirror 

Another illustration of the power of our microscopic 
model of optical elements is the parabolic mirror. In 
fact, it is possible to "build" out of two-level atoms mir- 
rors of arbitrary shapes. In the next simulation, the 
photon wave packet is propagating towards a parabolic 
mirror the shape of which is described by the equation 
X — xo + -^y^ = 2 — j^y^ ■ The focus of the parabola 
is at the point x — xq -\- ^ = —2.5, y = Q. The time 
evolution of the energy density is shown in Fig. pi At 
t = 0.0 the Gaussian photon is propagating towards the 
parabola. The little circle in between the photon and 
the parabola shows the position of the focus. At i = 8.0 
we see the photon wave packet being reflected from the 
parabola. We see the interference between the incoming 
and the re-emitted radiation. We note that at i = 12.8 
most of the radiation goes through the focus. In the last 
figure [t = 18.0) the photon wave packet propagates to 
the left. 
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FIG. 5. The energy density of the electromagnetic field at 
the moment when the incoming wave packet interfere with the 
radiation re-emitted by the atoms of the mirror (left) and the 
beam splitter (right). The central wavelength of the photon 
wave packet in the case of the mirror simulations is taken to 
be longer compared to the case of the beam splitter simula- 
tions. The interference pattern in the two cases is different. 
We see that in the case of the beam splitter part of the radi- 
ation is transmitted. The parameters of the simulations are 
specified in previous figures. 
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FIG. 6. The time evolution of the energy density of the 
one-photon wave packet reflected by a parabohc mirror com- 
posed of two-level atoms. At time i = the wave packet is 
localized to the left of the focus of the mirror (a little circle 
in the figure) and propagates towards it. At time t = 8.0 we 
see an interference pattern due to interference between the 
incoming wave packet and the re-emitted radiation. At time 
t = 12.8 the original wave packet is completely reflected by 
the mirror and is localized around the focus. The spatial de- 
pendence of the energy density is determined by the shape of 
the mirror. We can observe a reduction of the width of the 
reflected wave packet in the y direction. At time t — 18.0 
the wave packet is spread significantly. We see that the max- 
imal energy density is now smaller than in the original wave 
packet (compare with figure i — 0.0). The number of atoms 
from which the parabolic mirror is composed is Na ~ 1100. 
The parameters of the ato ms ar e the same as in Fig. H. The 
initial photon is Gaussian (3.12) with parameters a;o = —6.0, 
yo = 0.0, fc^o = 5.0, kyo = 0.0, AL = 0.125 and A|y = 0.125. 
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E. Interferometer 



FIG. 7. The time evolution of the energy density of the 
one-photon wave packet in an interferometer. The distance 
of the two mirrors from the beam splitter is the same. We 
see that the interference results in a wave packet propagat- 
ing upwards. The initial wave packet, the beam splitter and 
the mirrors have the same parameters as in figures considered 
above. Here the mirrors and the beam splitter are specified 
in Figs. and 0, respectively. 

On the other hand due to the quantum interference 
between the components of the radiation field coming 
from the two mirrors we can observe something com- 
pletely different: If the optical paths of the two compo- 
nents are equal then their relative relative phase is such 
that quantum interference results in an emergence of a 
single-photon wave packet traveling up (t = 45.0). On 
the contrary, if the distances of the two mirrors from the 
beam splitter are not equal then the relative phase of the 
two components which interfere on the beam splitter af- 
ter being reflected by the mirrors can result in a wave 
packet traveling left (see Fig. H). Here the difference of 
the optical paths is approximately one central wavelength 
of the original wave packet. We see that in this case most 
of the energy travels in a form of a wave packet to the 
left. 



Using a quantum beam splitter and two quantum mir- 
rors we can "construct" a single-photon interferometer 
(see Fig. 0). Here the one-photon wave packet comes to- 
wards the beam splitter (i — 0.0) and is divided into two 
parts which propagate towards the mirrors (t = 18.0). 
The distances of the mirrors from the beam splitter are 
exactly the same. The mirrors reflect the radiation back 
to the beam splitter. At t = 33.3 the two reflected parts 
reach the beam splitter. Each wave packet considered in- 
dividually would be splitted by the beam splitter into two 
parts going left and up (i.e. transmitted and reflected). 
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F. Two-slit experiment 
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FIG. 8. Same as in Fig. M except the distances of the mir- 
rors from the beam splitter differ by one wavelength. This 
difference leads to a quantum interference which results in a 
wave packet propagating to the left. The initial wave packet, 
the beam splitter and the mirrors have the same parameters 
as in figures considered above. 
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The microscopic quantum model wc study in this paper 
can be also used to study the two-slit experiment. Let us 
assume the photon wave packet which has a very broad 
energy density in the y-direction, i.e. this wave packet 
models a plane wave which approaches the mirror with 
two slits, see Fig. S. On the left we have placed another 
mirror. Without it, the part of the plane wave which is 
reflected from the double slit mirror would disappear at 
the left and reappear on the right because of the periodic 
boundary conditions we use in our simulations. 

The original one-photon wave packet (t = 0.0) prop- 
agates towards the mirror with two slits. At t — 5.0 
the "plane" wave packet is reflected from mirror. Some 
of the energy propagates through the slits (i.e. there is 
a nonzero probability that the original one-photon wave 
packet can be transmitted through the mirror via the 
slits). We see that through each slit a part of the energy 
propagates to the right - the interference between these 
components of the electromagnetic field are clearly seen 
(see t = 20.0 and t = 25.0). 

The "plane" wave packet which has been reflected by 
the mirror with two slits is then reflected by the left mir- 
ror and then again by the the double slit mirror. Here 
part of the energy "goes" through the slits again form- 
ing a second, more complex, interference pattern (see 
t = 20.0). This process of bouncing of the original wave 
packet between two mirror continues and each time a 
fraction of the energy passes through the slits. 

It is interesting to compare the interference structure 
with the theoretical prediction derived within classical 
optics. The formula can be calculated using Huygens' 
principle |I^. According to this, every point at the slit 
can be considered to be a source of secondary wavelets. 
The total intensity profile is a superposition of these 
wavelets. Let us first consider a text book treatment 
of a one slit mirror, Fig. p^. 



FIG. 9. The time evolution of the energy density of a pho- 
ton wave packet in the two-slit experiment. A "plane" wave 
packet propagates towards the mirror with two slits. The 
mirror is composed of two level atoms as in previous figures 
except in this case there are two slits now. Another quantum 
mirror is considered to be located to the left of the two-slit 
mirror. This configuration is chosen to make sure that none 
of the original energy passes to the right due to the periodic 
conditions we imposed on the Schrodinger equation. To make 
the figure more transparent we use a logarithmic scale for the 
energy density of the field. The number of modes in this sim- 
ulation is 512 X 512. The total number of atoms used in the 
mirrors is Na ~ 7872. 
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which has been created to the right of the two-sht mirror 
during the first reflection of the original wave packet: 



L/2 



/(</>) 



/(r, (j))dr. 



(4.4) 



To neglect the contribution of the second reflection we 
take the lower bound of the integral over the polar coor- 
dinate r to be r„iin = 10. The theoretical prediction ( |4.3| ) 
and the intensity derived from our simulations (4^) are 
shown in Fig. O. Both intensities are normalized in such 
a way that their maximum is equal to unity. Near 9 = 
the agreement between the two results is very good. For 
larger values of 9 there is a difference between the two 
lines which is understandable because we are comparing a 
classical result with a numerical simulation of a quantum 
model with realistic features such as the nonzero width of 
a mirror composed of two-level atoms or the wave packet 
which is not a plane wave, etc. Taking these differences 
into account it is surprising that the two pictures coincide 
so well. 



FIG. 10. We show one slit in the mirror. The width of the 
slit is d. The optical path difference of the radiation coming 
from two different points at the slit is Aa;. 

The plane wave packet of the frequency lu is coming 
from the left towards a slit of a width d. The field 
strength on the right coming from a specific point at 
the slit is proportional to the phase factor e*''^^"'^*) . The 
phase difference of the radiation coming from two differ- 
ent spatial points in the slit is Ax = ysui9, if the dis- 
tance from the mirror is long enough. According to the 
Huygens' principle the total radiation is a superposition 
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For two slits of width d and a separation a we have two 
integrals 
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which gives for the intensity 
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On the other hand, we can use results of our numeri- 
cal simulations and evaluate the intensity of the radiation 
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FIG. 11. We present the intensity of the radiation field at 
the region behind the two-slit mirror. The dashed line cor- 
responds to the intensity derived from a classical model [see 
Eq.(4.4)] while the solid line is obtained from our numerical 
simulations based on purely quantum description of the pro- 
cess. We see a very good agreement between the two results 
for small 6. 



G. Decay of a two level atom 

Till now we have considered in our simulations that 
the field has been initially excited and all the atoms were 
initially in the ground state. Obviously, our model can be 
also applied to a situation when one of the atoms is ex- 
cited and the field is initially in the vacuum state (i.e., we 
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still restrict ourselves to the one-excitation subspace of 
the total Hilbert space). In this section we briefly discuss 
the problem of a spontaneous decay of a two- level atom in 
a two-dimensional cavity. We consider the atomic tran- 
sition frequency to be lo — 15.0 and the dipole constant 
IS D = 0.05. The atom is situated at the origin [x — 0.0, 
y — 0.0) of the two-dimensional cavity. The number of 
the field modes is 256x256=65536. In Fig. |lj we present 
the natural logarithm of the excitation probability of the 
atom as a function of time. From here we can conclude 
that the decay of the atom is approximately exponential 
with a decay constant F ~ 0.14. In Fig. 13 we present 
probabilities of the excitation of the modes k^ {ky=0). 
Because the direction of the constant dipole vector of 
the atom is chosen to be in the z-direction, the ampli- 
tude profile is the same on any line which goes through 
the origin of the momentum space. As expected for times 
large enough the modes with \kx\ = 15 are dominantly 
excited. In fact the peaks are not exactly at the reso- 
nance frequency, there is a small shift which is identified 
to be a Lamb shift (from the figures we cannot see this 
but the shift can be determined from numerical values 
obtained in the simulation). 



Modeamplitudes 




FIG. 13. We present probabilities of the excitation of the 
modes with ky = 0. We see that for times large enough the 
modes with \kx\ — 15 are dominantly excited, i.e. the field 
mode at frequencies close to the resonant transition frequency 
of the atom are most excited. 



Decay of a two level atom 




FIG. 12. The exponential decay of an excited atom. The 
atomic transition frequency and dipole constants are ui = 15.0 
and D — 0.05, respectively. The atom is positioned in the 
center of the square cavity of the linear dimension L — lOvr. 
We consider 256 x 256 = 65536 modes of the electromagnetic 
field. The probability to find the atom in the excited state 
is plotted in the logarithmic scale - the exponential character 
of the decay is clearly seen. The corresponding decay rate is 
r = 0.14. 



We plot the energy density of the one-photon wave 
packet emitted by the decaying atom in Fig. n4 Be- 
cause of the rotational symmetry of the problem we plot 
just one "cut" (y = 0) in the energy density as a func- 
tion of X. The energy density is presented for two times 
t = 4.0 and t = 12.0. At both times there is a peak in the 
center where the atom is positioned. This means that at 
these two times the atom still emits the radiation (which 
is in agreement with the chosen decay rate F — 0.14). 
We turn our attention to the fact that at t = 4.0 the 
energy density is nonzero only for |a;| < 4.0. Analogously 
for the time t = 12.0 the energy density is nonzero only 
for |x| < 12.0. This reflects the fact that the causality is 
preserved in our simple quantum-mechanical treatment 
of the decay of the two-level atom in the cavity. Here we 
have presented just few features of the decay, the com- 
plete description of the process deserves more detailed 
discussion. For instance, one might be interested on how 
the decay depends on the mode spectra, the position of 
the atom, what are the values of the Lamb shift, how 
the decay depends on the frequency cut-off, etc. We will 
address these questions elsewhere. 
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Energy Density 



time = 4.0 
time = 12.0 




FIG. 14. We plot the energy density ol the one-photon 
wave packet emitted by the decaying atom. The parameters 
of the atom are the same as in Fig. hd 



V. CONCLUSION 

We have shown the results of many quantum mechan- 
ical simulations with two level atoms and a one photon 
wave packet inside a two dimensional cavity. The initial 
basis vectors are restricted to admit only one excitation. 
Because a rotating wave interaction between the radia- 
tion and the atoms is used, basis vectors with more ex- 
citations acquire no excitation. For these kinds of states 
the special numerical technique which utilizes FFT (Fast 
Fourier Transform) may be used. Using FFT the simula- 
tions become orders of magnitude faster, allowing more 
modes and atoms to be included. 

The atoms are at fixed positions and it is possible 
to build complicated structures with different kinds of 
atoms. Several layers of atoms which are on resonance 
with the incoming radiation form a quantum mechanical 
mirror if the density of the atoms is high enough. The 
mirror may have an arbitrary shape. In our simulations 
usual flat and parabolic mirror were used. One layer of 
detuned atoms forms a beam splitter. We have shown 
that, using mirrors and beam splitters, it is possible to 
build complicated optical networks. As an example the 
time evolution of a photon in an interferometer was stud- 
ied. 

Usually the optical components are taken to be classi- 
cal objects which give boundary conditions to the quan- 
tum mechanical time evolution or determine the modes 
used in a quantization. In our simulations the whole 
system including beam splitters and mirrors is in a well- 
defined quantum mechanical state. In addition to the 
simulations shown in this paper is is possible to build 
more complicated networks of beam splitters and mir- 
rors. One interesting possibility is to build cavities of 



arbitrary shape and study the time evolution of the pho- 
ton intensity inside the cavity. It is also possible to use 
moving atoms in the simulations allowing moving beam 
splitters and mirrors to be built. One extension of the 
current model would be to take basis states with more 
than one photon excitation into account. However, the 
number of basis states with a given excitation increases 
so rapidly that it is unlikely to be possible to use the 
methods of this paper for fields of higher intensity. Thus 
all the phenomena of nonlinear optics require novel com- 
putational approaches. 
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